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Abstract. We use numerical simulations of a bead-spring model chain to investigate 
the evolution of the conformation of long and flexible elastic fibers in a steady 
shear flow. In particular, for rather open initial configurations, and by varying a 
dimensionless elastic parameter, we identify two distinct conformational modes with 
different final size, shape, and orientation. Through further analysis we identify 
slipknots in the chain. Finally, we provide examples of initial configurations of an 
“open” trefoil knot that the flow unknots and then knots again, sometimes repeating 
several times. These changes in topology should be reflected in changes in bulk 
rheological and/or transport properties. 
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1. Introduction 

One area of complex fluids concerns the motion and topology of elastic filaments in fluid 
flows, which is inspired by both natural and industrial phenomena, such as flagellar 
motion m Ei E], fluid-structure interactions, e.g. [3 E], and polymer processing, e.g. 
O Hi E- m E]. Such long flexible filaments can have non-trivial dynamics, e.g. drift 
across streamlines can occur, e.g. [Zl nn H2], and topologies. For example, knots are 
found in Brownian systems such as bacterial DNA, protein structures, and polymer 
chains, and non-Brownian systems such as ordinary string, elastic fibers and chains of 
linked beads [13] [14] . In this paper we report the time-dependent shapes of flexible 
non-Brownian filaments in a steady shear flow and identify conditions that allow an 
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“unknotting-knotting” transition, where an unknotted filament is later observed to form 
an (open) knot as it rearranges continually in the flow. 

For rigid non-Brownian particles in shear flow, when the Reynolds number is low, a 
straight elongated particle has a Jeffery orbit, which is the basis for many studies, 
and helps organize those initial states that have orientation and rolling about the 
vorticity axis from those that primarily align with and tumble periodically about the flow 
direction. For flexible or deformable objects, characterizing the shape of the particle in 
flow, and the corresponding dynamics, is more challenging. For example, experimental 
measurements of fiber dynamics commonly use DNA as a model polymer [15], and, on 
a larger length scale where Brownian effects are less significant, linked colloidal chains 
are also used as experimental models of a worm-like chain [ I6j HZ] and the cellular 
flagellum [T8] . In all of these cases, dynamics in a single plane are measured, but it is 
difficult experimentally to simultaneously track dynamics in both the shear and vorticity 
directions [19] . 

On the other hand, numerical simulations make possible the visualization of 
multiple viewing planes and the detailed analyses of shape and topology of flexible 
filaments. Previous work, some of which includes Brownian effects and some of which 
does not, has identified and tracked the tumbling of flexible filaments (e.g. polymers) in 
detail, relating the tumbling period and shape to polymer length, flexibility and shear 
rate, as well as rates of rotation both in and out of plane [T21I201I2T1 22] [23]. However, due 
to computational constraints, these simulations of fibers are generally bead-chains with 
about N = 50 beads, with the longest ones containing N = 100 beads; an alternative 
approach to some of these questions is to use slender-body theory, e.g. [Ml 125] . 

In recent years, knots in long chains has been studied more widely, including the 
observation of thermally driven knotting and unknotting of a polymer chain [26], the 
influence of chain tension on the knot, e.g. [27] [2S], and the influence of electric fields on 
knots Also, the untying of a knot in an extensional flow has been studied recently 
[30], and has the spirit of the present paper where we ask about the effects of flow. 
In contrast, for non-Brownian systems other kinds of forcing, such as vibration, can 
produce knotting of chains of linked spheres naE]. We are not aware of any similar 
studies, experimental or numerical simulations, highlighting conditions that lead to long 
non-Brownian chains forming knots in flows. 

Thus, while significant steps have been taken towards understanding the detailed 
dynamics of flexible non-Brownian fibers, to understand topology is much more difficult: 
it requires three-dimensional detail and very flexible, long and thin fibers as, in theory 
(when there is no flow), the probability of knotting in a random walk increases 
exponentially with the length of the walk [31]. Moreover, a smart choice of initial 
conditions is needed and an operational definition is needed for a knot in an open fiber. 

In this paper, we focus on the dynamics and topology of long, flexible, non-Brownian 
chains of beads in a steady shear flow and use numerical simulations to investigate 
whether knots in the fiber occur as a result of the flow. We will find that unknotted 
fibers are capable of forming (open) knots, and in some cases we document a sequence of 
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unknotting-knotting transitions. The numerical solution of this kind of problem requires 
tracking the motion of N beads, where hydrodynamic interactions between the beads 
mean that the dynamics of the topology of the object involve 3 N coupled nonlinear 
ordinary differential equations. Not surprisingly, such dynamics should be expected to 
be chaotic though we have not attempted in this paper to link our study of complex 
shapes and the unknotting-knotting transition to the underlying chaotic dynamics. Our 
results highlight new aspects of the dynamics of flexible filaments in flow. 

2. Simulation methods 

2.1. A fiber as a string of spherical beads 

We consider a long non-Brownian flexible fiber in steady undisturbed shear flow 
Vo = 7 ze x of a fluid with viscosity g (figure CD). The vorticity vector is in the y- 
direction. The fiber deforms owing to elastic and bending forces and is modeled as a 
chain of N identical spherical beads with diameter d (figure UK) P2 .32]. The parameters 
of these forces are k, the spring constant normalized by irr/dfi, A, the bending stiffness 
normalized by nr/d^fi, and £q, the equilibrium distance between the closest bead centers 
measured in units of d. All distances are measured in units of d, and time in units of 
l/y. Thus, there are four dimensionless parameters that characterize the chain: k, A, 
N and £q. We work close to the inextensible limit with \i$ — 1| <C 1 and k 1 and 
focus on the dynamics for a wide range of values of A. 

We assume that the particle-scale Reynolds number is much smaller than unity, 
i.e. Re = <C 1, where p is the fluid density and the fluid flow satisfies the 

quasi-steady Stokes equations, with no-slip boundary conditions at the bead surfaces. 
The equations are solved using the multipole expansion method, which accounts for 
lubrication corrections to speed up the convergence of the friction matrix as the order of 
the multipole truncation is increased. This solution yields the velocities of the beads, as 
explained in [12]. The algorithm and its numerical implementation HYDROMULTIPOLE 
are described in [331IM] . Also, the application of these tools to model dynamics of flexible 
fibers is outlined in m The time-dependent translational velocities of the bead centers, 
evaluated by the hydromultipole numerical code, are used by an adaptive fourth- 
order Runge-Kutta procedure, from which the positions of every bead are updated. 

We focus on the dynamics of long fibers (N fs 100 — 150), which, in flow, tend to 
fluctuate from elongated shapes to compact shapes with some of the beads coming close 
together. Spurious overlaps of beads can occur, which can significantly influence the 
time scales of the dynamics of close particles |35| . Therefore, we analyzed different 
procedures that allow the dynamics to be continued once an overlap occurs and chose a 
procedure that caused the least number of overlaps. If the distance between two bead 
centers r < 1, then the lubrication correction for this pair is evaluated for the rescaled 
configuration 1 + e^, where ei u b = 10 -4 is the integration accuracy per step, and the 
dynamics of the original positions are continued. The details of this original simulation 
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Figure 1: Coordinate systems and notation for the simulations, (a) A ‘candy-cane’ 
shape, which is one of two typical initial configurations studied in this work. The 
expanded view labels beads i and i — 1. Also, is the position vector of each bead 
center and t t is the vector connecting two adjacent bead centers, (b) A configuration 
with shading to indicate the chain’s ellipsoid of revolution. The shear flow is shown in 
the background. The major axis of the chain’s configuration is shown with an arrow, 
while the red and blue coloring of beads denote two halves of the fiber for visual clarity. 
The orientation is measured by the angles 6 and q f>, such that 6 measures the inclination 
from the vorticity axis ( y ) and (j) measures the orientation in the xz-plane, measured 
from the 2 -axis. 


approach for handling close particle pairs are provided in the next section. 

In the simulations reported here, the initial geometry is chosen first from a random 
rotation of a ‘candy-cane’ shape with N = 152, i.e. a straight chain that bends into a 
semi-circle at one end; see figure [Hand Appendix A The candy-cane shape is convenient 
for starting with an elongated initial configuration that breaks the symmetries associated 
with shear flow. Also, we have performed calculations with a trefoil as an initial 
condition with N = 99 (see Section [4] and Appendix B ). For each simulation we fix 
£ 0 = 1.02 and k = 500, and vary A systematically. We evolve the shape in a steady 
shear flow using the HYDROMULTIPOLE algorithm 


2.2. A simulation procedure for compact structures 

As introduced in the previous section, and will be clear when the many simulations in 
this paper are shown, long flexible fibers can form compact structures with some of the 
beads coming very close to other beads. As a result, it sometimes happens that surfaces 
of pairs of beads spuriously overlap. This complication will be faced by all numerical 
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methods and for the long-time dynamics of long flexible fibers, which are the focus of 
this paper, such artifacts cannot be avoided. The physical problem is that the friction 
and mobility coefficients are not defined for the distance between the bead centers less 
than their diameter (equal to one in our units). In the literature, there exist different 
methods that allow the dynamics to be continued once an overlap takes place (e.g. see 
the brief review in (35)). In this sub-section we present a new approach for the problem 
of overlaps; the reader interested only in the results of the fiber dynamics, including the 
unknotting-knotting transition that we have identified, can skip this sub-section. 

Not surprisingly, it is known that numerical procedures for handling the interaction 
of close pairs can sometimes significantly influence time scales of close-particle dynamics. 
In particular, examples are known of systems of particles where spurious overlaps 
increase the period of oscillations, even by a factor of two [35]. Therefore, it is important 
to choose a numerical treatment of overlaps that minimizes the time-dependent number 
of overlapping pairs and also prevents the particle surfaces from remaining overlapped 
for a long time. 

We have observed that spurious overlaps of beads sometimes can significantly 
modify compact shapes of flexible fibers, lead to their orientation along the y (vorticity) 
axis, and change properties of the dynamics. In particular, a larger number of overlaps 
leads to a larger fiber curvature and a smaller radius of gyration. Therefore, to address 
the main questions of this paper we found it essential to use an optimal procedure once 
an overlap takes place. After analyzing and comparing a few such methods, we proposed 
a new one. Here we provide a brief description of this procedure (called BWE), which 
we have applied to continue evolution of the system of equations in case of spurious 
overlaps of beads. Our new approach gives a small number of overlapping pairs of the 
beads, which, in addition, usually do not remain overlapped for a long time. 

The procedure is based on modification of the lubrication correction in 

the case when the beads i and j overlap. To explain the idea, we first briefly explain 
the general idea how the lubrication correction is constructed [33 , 34]. Truncating the 
multipole expansion at the order L, we evaluate the multipole approximation for 
the friction matrix of the fV-bead system. This step can be accomplished even if the 
distance ry,- between the centers of beads i and j is close to but smaller than one (here, 
i, j = 1,..., N). We need to add a lubrication correction to speed up the convergence 
of the multipole expansion, + d ( r N> . This correction is the sum of all pairwise 

contributions, k ^ /, 

8 {2 \r kl ) = C {2 \r k i) - CfVfci), (1) 

where r k i is the relative position vector of particles k and l. Here, C i ' 2 \ r ki) and C^ 2 \ r ki) 
are, respectively, the two-body friction matrices - the exact one and its multipole 
approximation of the order L. 

Now, consider a pair of overlapping beads i and j, i.e. r tJ < 1, for which elements of 
the matrix are not defined - they diverge at the contact. We overcome this problem 
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and evaluate the dynamics from the unchanged positions of all of the particles. The 
multipole approximation for the friction matrix of the iV-bead system and the 
lubrication corrections for non-overlapping beads are not modified. 

The key idea of the method is to take e/ u & equal to the integration accuracy per 
step, which in this work is equal to 10” 4 . Owing to this choice, the time of a spurious 
overlap is reduced in comparison to smaller values of e^. Indeed, a typical overlap 
depth is of the order of the accuracy. Therefore, by assigning to the relative velocity of 
the beads the new, quite large value, which corresponds to the gap size between their 
surfaces, i.e. ei u b equal to a typical overlap depth, we obtain the required result that the 
overlapped beads can be disconnected easily, i.e. in relative terms they quickly separate 
their surfaces from each other. 


3. Results: Elongated initial configurations 


Upon starting the shear flow, the fiber, which was initially in the candy-cane 
configuration, described in Appendix A, aligns approximately with the flow direction 
and undergoes a ‘tumbling’ or ‘yo-yo’ motion (figure |2K,d; the chain halves are 
colored for clarity), similar to observations in studies of flexible polymers in shear flow 
PSIEEIEZIEE]. We identify two distinct behaviors for the dynamics, depending on the 
value of A. To quantify the dynamics, we report the time-dependent approximate shape 
and orientation of the chain. 


3.1. Typical simulations and geometric characterization 

First, for a given configuration, we calculate the moment of inertia tensor [39] 

1 N N 

= 2*5 £ £ b" - A [C - rf] (3) 

*= 1 j =1 

with the Cartesian coordinates n,m = 1,2,3, and the labels i,j of the vertices of each 
of the N beads. The lengths a, b, and c denote the three semi-axes of the ellipsoid 
of inertia associated with a configuration (figure [1 ]d). These lengths are related to the 
three eigenvalues X n of T nm , i.e. a = \/3Ai, b = a/ 3A 2 , and c = \/3A 3 . The largest 
eigenvalue a gives the primary axis of the ellipsoid, and the orientation angles of the 
associated eigenvector, 6 and 0, are measured from the vorticity and z axes, respectively, 
as shown in figure [ljn We use a, b , and c to visualize the fiber shape: if a ~ fe ~ c, the 
configuration is roughly spherical, whereas if a 6 ~ c, the shape is prolate ellipsoidal, 
and if a « 6 > c, the shape is oblate spheroidal. 

Also, we calculate the radius of gyration R g of each configuration, i.e. R 2 g = 

N 

A (r* — r cm ) , by comparing the coordinates of each bead r* to the center of mass r cm . 

i —1 
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(a) (b) (c) (d) 

Figure 2: Dynamics of chains for different values of A. (a) Time evolution of a chain, 
labeled (i) — (u), with A = l. The first three frames show a tumbling chain aligned with 
the flow, as the red and blue ‘halves’ have switched. At later times the shape is oriented 
along the vorticity (y) axis, (b) For A = 1, radius of gyration R g , axis lengths a, b, 
and c, and the orientation of the main axis as measured by 6 and </> each plotted versus 
time. Close pairs denote the number of beads that are within a distance d tJ /d — 1.05 
of another bead, not counting nearest neighbors in the chain. The dashed vertical lines 
denote times that correspond with the shapes — (c) The same measures plotted 

as function of time for A = 3. (d) Time evolution of a chain, panels (/)— (V), with A = 3. 


Since R 2 g = trace T nm then R 2 g = A 1 + A 2 + A 3 0 U- With these measures {R g , a, 6 , c, 9, (j)}, 
we track a given chain’s size, shape, and orientation as it evolves in the flow. 

For a given initial orientation, specified in Appendix A, chains with lower values of 
A, e.g. A < 2 , behave differently than chains with higher values, A > 3; qualitatively, 
the former form tight coils aligned with and rolling in the vorticity direction, while the 
latter forms more open tumbling configurations, by analogy to the results of Ref. m- 
We present detailed results from each category, starting from the same initial ‘candy- 
cane’ geometry. 

First, we focus on a simulation of a more flexible chain where A — 1; the evolution is 
shown in figure [ 2 a] and movie [ 2 a] and the dynamics are quantified in figure [ 2 b] Initially, 
the chain extends and orients in the flow direction, resulting in a large radius of gyration 
and highly skewed axes lengths (a S> b, c); the orientation angles of the main axis are 
6 fa 7 t /2 and (j) ~ 7 t / 2 , i.e. the flow direction. As time continues, the two ends fold 
towards the center, R g decreases significantly, while a decreases and becomes comparable 
to b and c. At still later times, the configuration remains relatively constant in the form 
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of a prolate ellipsoid (a b, c ) parallel to the vorticity axis; also, (j) systematically 

increases, which is consistent with a ‘rolling’ motion. 

We next consider a simulation of a stiffer chain, where A = 3; see figure Eh,d and 
movie Ml The time dynamics is different than the case of A = 1 , e.g. the collapse of the 
initially extended chain leads to a less compact shape that does not resemble a prolate 
spheroid. For most times, the orientation angles of the main axis remain close to the 
flow direction. 

Finally, we examine the number of close pairs for both A = 1 and A = 3. For each 
time-step, we check the distance between each pair of beads; we consider beads a close 
pair when dij, the distance between beads i and j, is less than some threshold, and the 
beads are not the closest neighbors, i.e. i^j±l. Here we set the threshold dij/d = 1.05, 
slightly larger than the dimensionless equilibrium inter-bead spacing £$ = 1.02. We 
notice that the number of close pairs in the more flexible fiber [A = 1) is around 4 times 
greater than the number in the less flexible fiber 1 figure l2bl l2cj) . These observations 
are consistent with A — 1 yielding a tightly packed coil and A — 3 resulting in a less 
compact shape. 

Both visually, as seen in figure Ek,d and movies [2k, d, and quantitatively, through 
figure [2}d,c, a change in A can result in considerably different dynamics. A closer 
inspection of the configurations shows that the chains may become self-entangled. For 
example, in panel (V) of figure [2dl the blue end has passed through a loop formed at 
the opposite red end. Visual inspection shows that this is a slipknot (see section [373]) . 
but such shapes raise the question of whether or not simple flows can cause chains to 
tie themselves into knots, and how to identify and classify such configurations. 

3.2. Knot classification 

For a given configuration at a given time, we then study the topology of a chain. A knot 
is a topological state of a closed curve, where no set of simple line deformations, e.g. 
‘pushing’ a strand, without passing one strand through another, can change its state to 
that of an open circle, or an ‘unknot’. To classify knots, polynomials known as knot 
invariants, which uniquely identify each curve to a knot type, can be calculated [42] . 
Here, we calculate the Alexander polynomial of each configuration (see Appendix C ). 

Any open curve, i.e. a curve with free ends, can always be ‘untied’ by line 
deformations and cannot be considered knotted in a rigorous sense. However, these 
‘open knots’ still have characteristics that can be classified, and affect the properties of 
a chain, such as its shape, size, and relaxation modes [39]; the classification of open knots 
remains a subject of study. In general, an open curve must have its ends joined together 
to form a closed curve, from which the topological state can be calculated. There are 
several well-studied closure schemes [43], but as these methods can result in different 
topological states for the same open configuration, there is inherently ambiguity. 

We use the stochastic closure scheme m- For each open configuration, we draw a 
sphere centered at the center of mass of the chain, with radius ten times the radius of 
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gyration. A point is randomly chosen on this sphere [32], and we connect by a straight 
line the two ends of the open curve to this point to create a closed curve that can be 
analyzed for its Alexander polynomial. Repeating this procedure with different random 
points on the sphere (here, 1000 randomly chosen points), we build a spectrum of knots; 
i.e. each configuration of the flexible chain can be described as a superposition state of 
various knot types. 

This approach effectively provides a statistical measure of the entanglement of each 
chain, which we report with a knotted fraction kf , i.e. the fraction of the stochastic 
closures which result in a knot of any type. For an open flexible fiber, “a knot” 
corresponds to kf close to 1. 

As might be expected from visual analysis of the configurations shown in figure [2] 
for the initial candy-cane configuration, no knots have been observed for the time of the 
simulation even though the fiber shapes were often very compact at long times. The 
calculated values of kf were always close to zero for all time steps in the simulation. 

Because the chains we study are highly flexible, it is also necessary to understand 
the topology of sections of the chain. Thus, for each subchain, the knotted fraction kf is 
calculated using the stochastic closure scheme and the resulting values are combined into 
a knot matrix K(i,j), where i, j = 1 • • • N. Each point (i,j) of the matrix corresponds 
to starting and ending beads of the subchain, i and j, respectively. For example, the 
origin, at (1,152), refers to a subchain starting at i — 1 and ending at j = 152, i.e. 
the entire chain. Another point on the matrix, such as (70,152), refers to a subchain 
starting at i = 70 and ending at j = 152, which is the latter half of the full chain. 
Therefore, every cell of the matrix K(i,j ) corresponds to a subchain of the original 
configuration. 

The knot matrix is symmetric, K(i,j ) = K(j,i), because the subchain that starts 
from i and ends at j is just the same subchain that starts at j and ends at %. Therefore, 
it is sufficient to show values of the knot matrix only for i < j. 

3.3. Slipknots 

In figure l2dTV) . we observe a configuration that has the blue end of the chain threaded 
through a loop formed from the red half of the chain. As, intuitively, a knot requires 
the pulling of an ‘end’ through a ‘loop’, the way we might tie a knot in a shoelace, 
these types of configurations bear further scrutiny. By closing the entire chain with the 
stochastic closure scheme, these configurations are shown to be unknotted - the fraction 
of nontrivial closures is low, kf close to 0. 

With this visual motivation, next we studied the results of simulations at different 
values of A, keeping an eye out for possible knots. For A = 1.4, the knot matrix K(i,j ) 
at an intermediate time when an entangled configuration occurs is plotted in figure [3] 
The white cells have zero knotted fraction, while the darker cells have a higher knotted 
fraction. We note that the knot matrix is not symmetric with respect to reflection 
in the line y = 152 — x, because, in general, (i,j) and (152 — j + 1,152 — i + 1) 
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Figure 3: Left: Example of a knot matrix for A = 1.4, N = 152, t = 682.7. The axes 
denote the first and the last bead number of the subchain. The cell is colored according 
to its knotted fraction: kf — 1 is black, and kf = 0 is white, as denoted in the inset 
scale. A wire-like scheme of the knotted subchain corresponding to (70,152) is shown. 
Right: The knot on the subchain (red) is only a slipknot in the entire chain, with the 
subchain pictured in red and the other sections highlighted in blue. 


are different subchains, and one of them can be knotted while the other one is not. 
For example, the subchains (70,152) and (1,83) are not the same and have different 
topological characteristics (black and white colors, respectively). 

In plots such as figure [3l dark regions can be identified as knot cores [46]. In these 
regions, the corresponding subchains have kf approaching unity and are considered 
knotted, even though kf is close to 0 for the chain as a whole. For example, as evident 
in figure [3] there is a darker block corresponding to a subchain that starts at the bead 
70 and ends at the bead 152. This knotted subchain is shown highlighted in red in the 
accompanying image in figure [3] In the right panel of figure [3l the entire chain is shown, 
with the knotted subchain highlighted in red, and the other part highlighted in blue. 
Thus, we have identified, qualitatively and quantitatively, the formation of a slipknot 
in a flexible chain subjected to a steady simple shear flow. 


4. Results: Trefoils in steady shear flow 


Now that we have explored one extended type of initial condition, we explore a more 
“knot-like” initial condition. Thus, we ask the question: can shear flow untie a trefoil? 
We set an initial condition in a trefoil (3i knot) shape (with one bead missing), see 


Appendix B, and vary the initial orientation of this shape relative to the shear flow as 


shown in Fig. [4] Then, we evolve the chain for different relative bending stiffness A. 
By using the stochastic closure scheme and counting invariants, we can calculate and 
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Figure 4: Initial trefoil shape and five distinct orientations. 


plot the fraction knotted kf at each time step, i.e. how many of the stochastic closures 
result in non-trivial knots. A knot corresponds to kf close to 1. 

We used five different orientations of this shape, shown in figure 0Ji-e. A trefoil with 
N = 99 is shown in figure 0k and denoted as Jo = a. Also, the orientations Iq — b and 
J 0 = c (see figure 0)3 and c) were obtained by rotating the configuration J 0 = a along 
the x-axis by 7t/2 and ir, respectively; the orientations Iq — d, and Iq = e (see figure HU 
and e) were obtained by rotating the configuration along the z-axis by 7t/2 and 3n/2, 
respectively. 

4-1. Dynamics: An unknotting-knotting transition 

We performed 29 simulations starting from a trefoil shape at different orientations 
Jo = a, b, c, d, e shown in figure 01 and several values of the bending stiffness A. For 
certain I 0 and A, we observed an unknotting-knotting transition, by which we mean 


Table 1: The existence of an unknotting-knotting transition as characterized by a trefoil 
unknotting then knotting ( 1 ) or not (0), as a function of the stiffness A and initial 
orientation Iq. The table summarizes the 29 different simulations that were performed. 
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Time =317 




Figure 5: The unknotting-knotting transition for an initial trefoil shape, A = 2.2, 
N — 99. The time evolution of the chain is shown, see also movie [5l Here kf denotes 
the knotted fraction based on the stochastic closure scheme described in the text. 



Figure 6: The unknotting-knotting transition for an initial trefoil shape, A— 1.4, N = 99. 
The time evolution of the chain is shown (see also movie |6|). Here kf denotes the knotted 
fraction based on the stochastic closure scheme described in the paper. 


that the trefoil unties and then ties again. In Table [U we specify values of the bending 
stiffness A and the initial orientation Jo for which such a transition was (or was not) 
observed. In most cases, after a transient period, the trefoil unties but, in some cases, at 
later times the chain knots again, as characterized by the knotted fraction, kf , decreasing 
to very small values and then rising up to almost one. 

We have found 13 examples of the unknotting-knotting transition: for fibers with 
A = 1.4, initially oriented as Jq = b, and for fibers with A = 2.1,2.2,2.3, initially 
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Figure 7: The unknotting-knotting transition for an initial trefoil shape, A =2.1, N = 99. 
The time evolution of the chain is shown (see also movie 0. Here kf denotes the knotted 
fraction based on the stochastic closure scheme described in the paper. 



Figure 8: The unknotting-knotting transition for an initial trefoil shape, A — 2.3, N = 99. 
The time evolution of the chain is shown (see also movie |HJ) - Here kf denotes the knotted 
fraction based on the stochastic closure scheme described in the paper. 


oriented as J 0 = c. In most cases, the unknotting-knotting transition is observed more 
than once during evolution of the same fiber. In figures [MSI we plot the time-dependent 
knotted fraction kf based on the stochastic closure scheme, and display some unknotted 
and some knotted configurations. Also, in the figures we use different colors to help 
highlight the complex topology of the complicated chain. 

Two examples of the unknotting-knotting transition with A = 2.2, visible in the 
time evolution of kf in figure |5l are shown in movie [5] In particular, in figure [5] it is also 
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illustrated how the chain first unknots as kf decreases below 0.1, and then knots again, 
as evidenced by kf then increasing to approach unity. Eventually, however, the chain 
is visibly unknotted and kf remains near zero. The other examples of the unknotting- 
knotting transitions that we have uncovered are shown in movies [6]|8] and the plots of 
the time-dependent kf in figures EHE1 Specifically, for A = 1.4, as displayed in figure [6] 
and movie [61 a knotted fiber is unknotted after time t ps 330 and knots again at time 
t fs 380. At the end of the simulation, the chain remains in the knotted configuration, 
with kf close to one. The evolution of a fiber with A = 2.1 is shown in figureOand movie 
[0 Again, the initially knotted chain unknots and knots again several times. Finally, 
the chain is knotted and stays in this form to the end of the simulation. 

In yet another example, the dynamics of a fiber with A = 2.3 is presented in the 
figure [8] and movie 0 The fiber unties and knots again many times, as shown in the 
graph of the knotting fraction kf. Finally, the fiber unties and stays unknotted until the 
end of the computation. Thus, we have examples, at least for the time duration of our 
simulations, of some cases that end up knotted and other cases that end up unknotted. 

All of the configurations displayed in figures [MHl can be seen in three dimensions, 
and rotated using the corresponding MATLAB fig-files, given in the Supplementary 
Data and labeled to match the numbers of the related figures; for example, the hie 
hg5a.£g corresponds to the first frame of figure 5 in the text. 


4-2. Characterizing configurations that display the unknotting-knotting transition 

Comparing the different simulations with each other, e.g. figures [MHl we conclude that 
the existence of the unknotting-knotting transition is not correlated with a final knotted 
or unknotted state. The question then can be asked if the existence of the unknotting- 
knotting transition for a trefoil is correlated with some specific mode of the flexible 
fiber dynamics, similar to those observed in the simulations that started from the initial 
candy-cane configuration, which were shown in figured For example, the characteristic 
features of the two main conformational modes where the fiber primary axis tends to 
be aligned along the vorticity axes (mode I) or the how direction (mode 11), discussed 
earlier in connection to figure [2j are summarized in table [2l 


Table 2: Characteristics of two main shape modes of hexible hbers in shear how. 


parameter\mode 
a, b, c 

0 

<f> 

close pairs 


I 

a » b > c 
ps 0 or 7r 
increasing 
many 


II 

a>6»c 

PS 7t/2 

■- ZL nr ( 2 k~H-) 7r 
' 2 2 

a few 


Therefore, in hgure 0 we compare the dynamics of trefoils with different bending 
stiffnesses A and initial orientations J 0 in terms of their characteristic modes of 
conformation. We present evolution of the same parameters of the hber shapes as 
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Figure 9: Dynamics of trefoils for different values of bending stiffness A and initial 
orientations Jo- The fiber radius of gyration R g , length of the semi-axes a, 6, c, 
orientation 6 and </> of the primary axes, and the number of close pairs are defined in 
the main text. The unknotting-knotting transition is observed during evolution shown 
in columns (ii) and (iii). 


in figure El but this time for fibers that have an initial trefoil shape. For two sets 
of the parameters, corresponding to the dynamics shown in columns (ii) and (iii), the 
unknotting-knotting transition is observed. For two other sets of the parameters, related 
to columns (i) and (iv), the unknotting-knotting transition is not observed. 

We first compare two cases for which the unknotting-knotting transition is observed, 
i.e. columns (ii) and (iii) in figure |9] For example, the evolution displayed in column 
(iii) for larger times approximately corresponds to mode I, which is characterized by 
an elongated shape with relatively large aspect ratios, and the tendency to align with 
the vorticity direction and roll around it. For comparison, at long times, the dynamics 
visible in column (ii) approximately corresponds to mode II, which is characterized by 
a smaller elongation, and no special features of the dynamics related to the vorticity 
direction. 

On the other hand, comparing columns (i) and (ii) in figure [9l we conclude that 
the same mode II configuration is reached by fibers for which the unknotting-knotting 
transition, respectively, did not and did take place. These examples of two topologically 
different behaviors correspond to different initial orientations of the trefoil characterized 
by a single value of A. The complementary comparison of columns (iii) and (iv) in figure 
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0 demonstrates that fibers with and without the unknotting-knotting transition can both 
later exhibit the same mode I configuration. In this case, examples of two topologically 
different behaviors correspond to the same initial orientation of the trefoil but different 
values of A. Therefore, the examples displayed in figure [9] illustrate that the existence of 
the unknotting-knotting transition is not correlated with the final configuration mode I 
or mode II of the dynamics. 

Finally, we investigate if the process of creating a knot is correlated with any 
simultaneous specific features of the fiber shape. Two examples shown in figure [TO] 
contribute to a negative answer. In both cases, the fiber starts from very small knotted 
fraction kf, which increases up to almost one. However, in the first case (shown in the 
left part of figure ITOT) . the fiber is in mode II, and in the second case (shown in the right 
part of figure flOl) . the fiber is in mode I. To conclude this discussion, we note that the 
process of knot creation takes place both for fiber shapes that are compact (figures [6] 
[7]) and loose (figures [5] |HJ) • 


,4 = 1.4 Jo = 6 4 = 2.1 J 0 = c 



Figure 10: The process of creating a knot can take place when the fiber is in mode I 
(right) or II (left). 
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We have focused on the question of whether a simple shear flow can tie a knot in a 
flexible fiber. If yes, then can we predict the future existence of a knot based on the 
fiber flexibility and its initial shape and orientation? Also, is the knot a transient 
effect or a stable configuration? To address these questions, the fiber was modeled as 
a chain of N beads connected by very stiff springs with additional elastic forces acting 
against the fiber bending. The fiber evolution was determined using the numerical 
code HYDROMULTIPOLE, based on the accurate multipole algorithm for solving the 
Stokes equations with no-slip boundary conditions on the bead surfaces. In addition, to 
characterize the topology a knot on a flexible open fiber was defined as in Refs. [33-34] 
and associated with the knotted fraction kf close to unity. 

In the first part of our study, we investigated the evolution of fibers of different 
flexibilities, but with the same initial orientation and the same untied shape (with a 
vanishing knotted fraction). We concluded that for such a simple initial configuration, 
a change of flexibility is not sufficient to cause formation of a knot, although evidence 
of slipknots was observed. 

Then, in the second part of our study, we constructed a more complicated initial 
shape: a trefoil with the knotted fraction equal to one. We used the same initial shape 
at different initial orientations and for fibers of different flexibilities, and we followed 
their evolution in a shear flow. A few of the more flexible fibers remained knotted, but 
most of the fibers untied after some time. In particular, we then observed that some 
of these untied fibers again became knotted, as illustrated in figure fill for A = 2.3 and 
N = 99. In this figure we display two shapes with the knotted fraction close to zero, 
which later on lead to evident knots (with the knotted fraction close to one). 

Therefore, indeed, a simple shear flow can tie a knot on a flexible fiber. Sometimes, 
such an unknotting-knotting transition repeated several times during the evolution of 
the fiber dynamics. We observed 13 such unknotting-knottiug transitions in 4 (out of 
29) simulations, for 4 different (moderate) flexibilities and 2 different initial orientations. 
Based on all of our observations, we conjecture that the formation of knots is a transient 
effect and may possibly be related to the chaotic properties of the dynamics of this N 
bead system. 

To conclude, we have analyzed basic features of the dynamics of long flexible fibers 
in unbounded steady shear flow. In particular, we give examples of flexible fibers, ini¬ 
tially shaped as open trefoils, which are later untied and then knotted again, sometimes 
several times, by the flow. This finding provides a new perspective for future studies 
of topological structures of long flexible non-Brownian chains evolving as they are en¬ 
trained by ambient flows, and possibly changing the macroscopic transport properties. 
Finally, as we noted in the introduction, the numerical solution of this problem requires 
the solution of 3N coupled nonlinear ordinary differential equations for the positions of 
the N beads. We should expect such a dynamical system to have all of the attributes 
of chaos. How the observed dynamics of knot formation and unknotting may (or may 
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Figure 11: From unknotted to knotted: examples that illustrate that simple shear flow 
can tie a knot in a flexible fiber. Knotted fraction close to zero (left) is followed by 
knotted fraction close to one (right). A = 2.3 and N = 99. 


not) link to the underlying chaotic trajectories is a subject for future investigations. 
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Appendix A. Initial configurations of a “candy-cane” shape 

In the first part of our study, we started simulations from initial configuration of a fiber 
made of 152 beads, shown in figure IA11 This ‘candy-cane’ shape was generated from 
the equation (20 coss, 20 sins + 20, 0) with s ranging from —7 t/ 2 to 7 t/ 2 and the step 
ds = 0.05105. This equation is a semicircle of radius 20, and ds has been chosen in 
such a way that the distance between each point is equal to £ 0 — 1-02 bead diameters. 
A rod starting from (—91.8,0,0), ending at (0,0,0) is then connected; the points are 
again £ 0 = 1-02 apart, and described by: 

[-1.02(91-s), 0,0], (A.l) 

with s ranging from 1 to 90 in increments of 1. A rotation matrix is randomly generated 
and applied to these coordinates to rotate them out of the plane xy. In figure IA11 we 
show the initial condition which is used as an input for the simulations displayed in 
figure 2 and discussed in this paper. 
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Figure Al: Initial candy-cane shape. 


Appendix B. Initial configurations of a trefoil shape 

The trefoil shape is constructed with the use of the following equation, 

[sin s + 2 sin(2s), coss — 2 cos(2s), — sin(3s)], (B.l) 

from which 100 evenly spaced points are chosen with i 0 = 1.02. Since the points are 
evenly spaced along the path of the curve, there are slight variations in the straight-line 
distance between each point. One point was then removed to ‘open’ the trefoil as an 
open curve. The resulting trefoil with N = 99 and kf ~ 0.98 is shown in figure [4^, and 
denoted as Jo = a. 

Appendix C. Knot Invariants 

After our open curve is closed by the stochastic closure scheme, it can be analyzed 
for its polynomial knot invariant. There are by now many different knot invariants, 
which have varying ability to distinguish knots; the first one, which was discovered in 
1928 by Alexander, is one of the simplest to compute. It has some drawbacks, such as 
the inability to distinguish handedness (i.e. it cannot distinguish between two mirror 
images), but as we only need to distinguish between an ‘unknot’ and a ‘knotted’ curve 
to calculate the knot fraction, the Alexander polynomial is sufficient. 

The invariant is calculated by projecting a curve into 2D, then identifying ‘crossing 
points’ and ‘regions’, while keeping track of the ‘overlying’ and ‘underlying’ branches at 
each crossing. If there are n crossings, there will be n + 2 regions in the knot diagram. 
Depending on the orientation of the overlying/underlying branches at each crossing, the 
regions adjacent to the crossing are assigned different indices. A n x (n + 2) incidence 
matrix is constructed from the indices at each crossing and region, and two columns 
corresponding to adjacent regions are removed, to form a nxn matrix. The determinant 
of the matrix is then the Alexander polynomial. More details, and proof that removing 
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two arbitrary columns (provided they refer to adjacent regions) results in the same 
normalized Alexander polynomial, can be found in his 1928 paper m 

We have written a custom code to compute the Alexander polynomial of any 
configuration, given the (x, y, z ) coordinates of each bead in the chain, following the 
algorithm described by Alexander, with the help of functions from other Alexander 
polynomial calculators [l8j . 
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